A new non-perturbative method to compute accurate os- 
cillation modes in rapidly rotating stars is presented. In this 
paper, the effect of the centrifugal force is fully taken into ac- 
count while the Coriolis force is neglected. This assumption is 
valid when the time scale of the oscillation is much shorter than 
the inverse of the rotation rate and is expected to be suitable 
for high radial order p-modes of S Scuti stars. Axisymmetric 
p-modes have been computed in uniformly rotating polytropic 
models of stars. In the frequency and rotation range considered, 
we found that as rotation increases (i) the asymptotic structure 
of the non-rotating frequency spectrum is first destroyed then 
replaced by a new form of organization (ii) the mode ampli- 
tude tends to concentrate near the equator (iii) differences with 
\^ perturbative methods become significant as soon as the rota- 
^D tion rate exceeds about fifteen percent of the Keplerian limit. 
The implications for the seismology of rapidly rotating stars 
are then discussed. 
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Abstract. 



1. Introduction 

Since helioseismology revolutionized our knowledge of the so- 
lar interior, great advances in stellar structure and evolution 
theory are expected from asteroseismology. Major efforts in- 
cluding space missions are under way to detect pulsation fre- 
quencies with unprecedented accuracy across the HR diagram 
(Catala et al., 1995; Walker et al., 2003). To draw information 
from the observed frequencies, seismology relies on the theo- 
retical computation of eigenmodes for a given model of star. 
Yet, except for slowly rotating stars, the effect of rotation on 
the gravito-acoustic modes is not fully taken into account in 
the present theoretical calculations (e.g. Rieutord 2001). 

Rotational effects have been mostly studied through per- 
turbative methods. In this framework, both O/w, the ratio of 
the rotation rate Q to the mode frequency o), and Q/ -\JGM/R^, 
the square root of the ratio of the centrifugal force to the grav- 
ity at equator are assumed to be small and of the same order 
Solutions valid up to the first, second, and even third order 
in Q / ai hav e be en obtained respectively by Ledoux (1951), 
ISaiol ( Il98lh and ISoufi et alJ (|l998). The first order analysis 
proved fully adequate to match the observed acoustic frequency 
of the slowly rotating sun (Dziembowski & Goode, 1992). At 
the other extreme, the perturbative methods are not expected 
to be correct for stars approaching the Keplerian limit Q.]^ = 
■\JGM/Rl, where Re is the equatorial radius. Achernar is a spec- 
tacular example of such star since interferometric observations 
showed a very important distortion of its surface, the equato- 
rial radius Rg being at least one and a half time larger than the 
polar radius Rp (Domi ciano de Souza et all 120031) . In the con- 
text of Roche models, such a flattening occurs at the Keplerian 
limit Q.X- For intermediate rotation rates, second or third or- 
der perturbative methods might be used, but the main prob- 
lem is that the limit of validity of the perturbative methods is 
unknown. Departures from the perturbative results would im- 
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pact the values of individual frequency but also other proper- 
ties which are commonly used to analyze the spectrum of ob- 
served frequencies. This concerns in particular the rotational 
splitting, the asymptotic large and small frequency separations 
or the mode visibility. 

New methods able to compute accurate eigenmodes in ro- 
tating stars are therefore needed to allow progress in the seis- 
mology of rapidly rotating stars. Incidentally, such methods 
would also assess the limit of validity of perturbative analy- 
sis. The main difficulty comes from the fact that, except in the 
special cases of spherically symmetric media and uniform den- 
sity ellipsoids, the eigenvalue problem of gravito-acoustic res- 
onances in arbitrary axially symmetric cavities is not separable 
in the radial and meridional variables. For self-gravitating and 
rotating stars, a two-dimensional eigenvalue problem has to be 
solved. 

IClemenll ( ll98lLll998l) made the first attempts to solve this 
eigenvalue problem for gravito-acoustic modes, investigating 
various numerical schemes. However, the accuracy of his cal- 
culations is generally difficult to estimate. Moreover, the differ- 
ent numerical schemes could not converge for low frequency g- 
modes when Q/w exceeds about 0.5. Since then, eigenmodes 
in this frequency range have been successfully calculated by 
Dintrans & Rieutord (2000) using spectral methods. These au- 
thors however did not consider the effect of the centrifugal ac- 
celeration in their model. The search for unstable modes in neu- 
tron stars also triggered the development of numerical schemes 
able to solve the two-dimensional eigenvalue problem. But 
only surface gravity modes (f-modes) and some inertial modes 
(r-modes) have been determi ned in this context ( Yoshida et al.l 
l2005l) . lEspinosa et alJ ( 120041) reported calculations of adiabatic 
acoustic modes in MacLaurin spheroids of uniform density ne- 
glecting the Coriolis force, the potential perturbation and the 
Brunt- Vaisala frequency. 

In this paper, we present a new method to compute ac- 
curate eigenmodes in rotating stars. For the first application 



F. Lignieres et al.: Acoustic oscillations in rapidly rotating polytropic stars 



of the method, we only consider the effect of the centrifugal 
force through its impact on the equilibrium state of the star; 
we thus neglect the Coriolis force. This assumption is valid 
when the time scale of the oscillation is much shorter than the 
inverse of the rotation rate and is expected to be suitable for 
high radial order p-modes of 6 Scuti stars. The problem is fur- 
ther simplified by using uniformly rotating polytropes as equi- 
librium models and assuming adiabatic perturbations as well 
as the Cowling approximation. Low frequency axisymmetric 
p-modes have been computed for rotation rates varying from 
Q = up to Q/Q/f = 0.59, this ratio corresponding to a typical 
6 Scuti star (M = 1 .8Mq, R - 2Rq) with an equatorial velocity 
of 240 kms"' . The centrifugal force modifies the effective grav- 
ity in two ways: it makes it smaller and aspherical. Decreasing 
the eff'ective gravity should affect sound waves by reducing the 
sound speed inside the star and by increasing the star's volume, 
thus potentially the volume of the resonant cavity. Besides, the 
physical consequences of the non-spherical geometry are basi- 
cally unknown. 

In the following, the formalism and the numerical method 
are presented along with accuracy tests. Then, the parameter 
range of the calculations is given together with the method used 
to label the eigenmodes. The structure of the frequency spec- 
trum, some properties of the eigenfunctions and the differences 
with perturbative methods are further analyzed as a function of 
the rotation rate. These results are discussed in the last section. 

2. Formalism 

Accurate numerical solutions of 2D eigenvalue problems re- 
quire a careful choice of the numerical method and the math- 
ematical formalism. In this section we explain the choices that 
have been made for the variables, the coordinate system, the 
numerical discretization, and the method to solve the resulting 
algebraic eigenvalue problem. All play a role in the accuracy 
of the eigenfrequency determinations that will be presented at 
the end of this section. 



2.1. Equilibrium model 

The equilibrium model is a self-gravitating uniformly rotating 
polytrope. It is therefore governed by a polytropic relation, the 
hydrostatic equilibrium in a rotating frame, and Poisson's equa- 
tion for the gravitational potential: 

l + l/N 



= -VPo-poV (^0-^2*2/2) 



(1) 

(2) 
(3) 



where Pq is the pressure, po the density, K the polytropic con- 
stant, A^ the polytropic index, i^o the gravitational potential, s 
the distance to the rotation axis and G the gravitational con- 
stant. 

The polytropic relation and uniform rotation ensure that the 
fluid is barotropic. A pseudo-enthalpy can then be introduced 
ho - J dPo/po = (1 + N)Po/po and the integration of the hy- 
drostatic equation reads: 

ho^hc- (i//o - ijfc) + -^^s^ (4) 



where the subscript "c" denotes the value in the center of the 
polytropic model. Equation is then inserted into Poisson's 
equation to yield: 

A^„=4.Gp.(l-^:^.^f (5) 



he 2hc I 

Equation (|5} is solved numerically, using an iterative 
scheme. Since the shape of th e star is not spherical, a s ystem of 
coordinates ({,6,(f>) based on Bon azzola et alJ (11998) is used, 
such that ^ = 1 corresponds to the surface of the spheroid (more 
details on the coordinate system are given in section 2.3). This 
enables the use of spectral methods both for the radial coordi- 
nate ^ and the angular ones. The numerical method is further 
detailed in Rieutord et al. (2005). 

2.2. Perturbation equations and boundary conditions 

Neglecting the Coriolis force, the linear equations governing 
the evolution of small amplitude adiabatic perturbations read: 



d,p + V- (pov) = 0, 

Pod,v = -VP + pgo - poViff, 

d,P + v-VPo^cl (5,p + v-Vpo), 

Ai^ - AnGp 



(6) 
(7) 
(8) 
(9) 



where v, p, P and \p, are the perturbations of velocity, den- 
sity, pressure and gravitational potential. The sound speed is 
Co - yj^i.aPo/po, Ti^o being the first adiabatic exponent of the 
gas, and the effective gravity go - -V (i^o - Q.^s^/2] has been 
introduced. In the framework of polytropic models of stars, the 
polytropic relation Q is assumed to give a reasonably good ap- 
proximation of the relation between the pressure and the den- 
sity of the equilibrium state. As the first adiabatic exponent Fij) 
is obtained from the equation of state of the gas, Fio is in gen- 
eral not equal to 1 -H 1 /N. 

We simplified Eqs. (|6}, Q, (JHll, © following two con- 
straints: First, the governing equations should be written for 
general coordinate systems because we shall use a surface- 
fitting non-orthogonal coordinate system. Second, they should 
take the form A1Z - AQX where A is the eigenvalue, X is 
the eigenfunction, and Al and Q are linear differential opera- 
tors. Indeed, the method that we shall use to solve the alge- 
braic eigenvalue problem obtained after discretization works 
for problem reading [M]X - A[Q]X, where X is the discretized 
eigenvector and, [M] and [Q] are matrices. Taking the time 
derivative of Eqs. Q and (|5) and using Eqs. (|6j and (|SJ to 
eliminate the pressure and density perturbations, we obtain two 
equations for the velocity and the gravitational potential pertur- 
bation: 



dlv = V (c^ -H V ■ go - d,iff) -xAo 
Ad,ip - -AnG (v ■ Vpo + pQx) 



(10) 
(11) 



where ;if = V ■ v is the divergence of the velocity. The vector Aq 
characterizes the stratification of the equilibrium model: 



Ao 



1 Vfo 

r.o Po 



Vpo 
po 



fo^o'. 



r«o, 



(12) 



F. Lignieres et al.: Acoustic oscillations in rapidly rotating poly tropic stars 



where A^o is the Brunt-Vaisala frequency and mq is the unit vec- 
tor in the direction opposite to the effective gravity defined by: 



^0 



-||^oll«o. 



(13) 



Note that the barotropicity of the fluid has been used to obtain 

In addition to the gravitational potential perturbation, the 
right hand sides of Eqs. ( I10> and (II 1> only involve the diver- 
gence of the velocity and the scalar product of the velocity with 
vectors parallel to gravity. Then, the scalar product of Eq. ( I10> 
with gravity. 



dlv ■ go = go ■ V (cIy + v- go- dtijj) -xgQ ■ ^d 

and the divergence of Eq. ( llOl l. 

dlx = A(cf^ + V ■ go - d,t//) - V ■ OfAo) 



(14) 



(15) 



yield, together with Eq. Jlll l. a complete set of differential 
equations for the variables v ■ go, X ^"d tfr. Pekeris ( 1938) who 
studied the oscillations of spherically symmetric polytropes 
considered similar variables but, instead of Eq. ( I10> . used a 
combination of Eqs. (I14t and (I10> which does not involve sec- 
ond order derivative with respect to the radial coordinate. For 
general system of coordinate as well, the order of the differen- 
tial system can be lowered replacing Eq. JlOt by the following 
one: 



X-8''di 



v-go 
8o )l 



£(clx + V ■ go - d,>J/) ■ 



V-(xAo)+g''di 



Xgo ■ Ao 



where the linear operator X,, defined by, 



£(F) = AF-g''di 



go-VF 



8i 



(16) 



(17) 



does not contain second order derivatives with respect to the 
first coordinate x'. In this expression, g^ is the first contravari- 
ant component of the gravity in the natural basis {E\,E2,E-i) 
defined by Ei - dOM/dx', and g" is the first contravariant 
component of the metric tensor. 

The equations are non-dimensionalized using the equatorial 
radius, Rg, as length unit, the density at the center of the poly- 
trope, pc, as density unit and To - (47rGpc) '^^ as time unit. As 
we seek for harmonic solutions in time, the variable are written 
F - Fexp(ia)t). Dropping the hat and denoting dimensionless 
quantities as previous dimensional ones, the governing equa- 
tions read: 



AW = go ■ V (c^ + W + ^) - clN^x 

£(clx+W + ^')-V-(xAo)+g''di 
= AH' - doW - pox 



^_Xx^ 
8o ) 



(18) 

(19) 
(20) 



where A - ■ 
IIVpoll 



-w , W — V ■ go, ^ - -iojip and do denotes 



do = 



llgoll 



(21) 



Another form of these equations may be obtained replacing W 
by a new variable U - c^ + W + ^. The set of equations then 
reads: 



go-VU-clNlx^MU-^-clx) 



(22) 



j:{u)-v-{xAo) + g''di 



-^"5,( 



f»' 
4 



fZ-Tl 



' 8o 



c^ 



(23) 



^0 

- doU + (docl - po)x + fl'o*!' + A^ = (24) 

As in lPekerisI ( Il938h . the boundary conditions are that the grav- 
itational potential vanishes at infinity and that U and x be reg- 
ular everywhere. 

2.3. Coordinates choice 

The choice of the coordinate system has been guided by two 
considerations. First, for the accuracy of the numerical method, 
it seems preferable to apply the boundary conditions on a sur- 
face of coordinate. This imposes that the stellar surface is de- 
scribed by an equation ^ - constant, where f is one of the co- 
ordinates. Second, when using spherical harmonic expansions, 
the regularity conditions at the center have a simple form for 
spherical coordinates only. Therefore, the coordinate system 
should become spherical near the center If (r, 6, (p) denotes the 
usual spherical coordinates and r - S{6) describes the surface, 
families of coordina tes (4", 9' , </>') ve rifying both conditions have 
been proposed bv Bonazzola et al.l (il998) : 



K^,0') 



(25) 



where 



ri(,0) = a( + AiO[S(e)-a] (26) 

The polynomial A{() = (5(^ - 3^^)/2 ensures that r cc ^ near 
the center and that ^ = 1 describes the surface r - S(ff). The 
scalar a is chosen so that the transformation (r, 0, (p) i-> ((, 0, 0) 
is not singular and the numerical convergence is fast. We find 
that a = 1 - e is a convenient choice, e - I - Rp/Re being the 
flatness of the star surface. In the following, we shall refer to ^ 
as the pseudo-radial coordinate. 

To express the governing equations in this non-orthogonal 
coordinate system, we use the covariant and contravariant com- 
ponents of the corresponding metric tensor The non-vanishing 
components read: 

^? 

"^^ " ',2^ ^ ..2„ . (27) 



ni 



8i2 = ^21 = r^rg 
822 ^ r + r^ g33 = r^ sin^fl 

g'' = ('-' + rl)/(r^rp g'^ = g^i ^ ^rg/ir^r^) 



,22 



= 1/r^ 



g^' = l/(r^sm^e). 
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and the square-root of the absolute value of the metric tensor 
determinant is: 



'r^ sin0. 



(28) 



In Appendix A, the linear operators involved in Eqs. M2\ . 
\2'i\ . (I24> are expressed in terms of the partial derivatives of 
r(^, ff). Note that, for vectorial operators, we used the natural 
basis defined above. 

2.4. The numerical method 

The method follows and g eneralizes the one presented in 
iRieutord & Valdettarol (119971) . The numerical discretization is 
done with spectral methods, spherical harmonics for the an- 
gular coordinates 6 and cp and Chebyshev polynomials for the 
pseudo-radial coordinate (. The variables U, *P and x are ex- 
panded into spherical harmonics: 



e=Q m=-e 



(29) 



and equivalent expressions for ^' and x, where £ and m are re- 
spectively the degree and the azimuthal number of the spherical 
harmonic Y"'{0, (p). Then, the governing equations are projected 
onto spherical harmonics to obtain a system of ordinary differ- 
ential equations (ODE) of the variable ^ for the coefficients 
of the spherical harmonic expansion 11^1^(0, xi,(0' ^LiO- This 
system is then discretized on the collocation points of a Gauss- 
Lobatto grid associated with Chebyshev polynomials. It results 
in an algebraic eigenvalue problem [M]X = /l[2]^, where X 
is the eigenvector of Lx Nr components and [M] and [Q] are 
square matrices of L x A^, lines, L and Nr being respectively the 
truncation orders on the spherical harmonics and Chebyshev 
basis. The algebraic eigenvalue problem is solved using either 
a QZ algorithm or an Arnoldi-Chebyshev algorithm. The QZ 
algorithm provides the whole spectrum of eigenvalues while 
the iterative calculation based on the Arnoldi-Chebyshev algo- 
rithm computes a few eigenvalues around a given value of the 
frequency. 

Because of the symmetries of the equilibrium model with 
respect to the rotation axis and the equator, one obtains a sep- 
arated eigenvalue problem for each absolute value of the az- 
imuthal number | m \ and each parity with respect to the equator 
Thus, for a given m > 0, we have two independent sets of ODE 
coupling the coefficients of the spherical harmonic expansion 
having respectively even and odd degree numbers, that is: 






(30) 



where Q < k < -t-oo. The solutions of these two ODE systems 
are either symmetric or antisymmetric with respect to the equa- 
tor 

The two sets of ODE can be written in the form: 



^d^ 
^W"'- 



J— -HCId S 



d( 



£ S, 



(31) 



where E denotes 
4,+ 



or a 



u- 

X' 

4*- 



(32) 



and where the matrices are defined by blocks as follows: 



y[: 



c 



£) = 



CO 


\ 








I A33 j 



fBii 
B21 B22 
B33; 



rcii C12 

C21 C22 
C31 C32 C33 ^ 



i 

D21 D22 







-D21 





£ = 



-11 E12 


-Ell 


12 1 E22 


-E21 









(33) 



(34) 



(35) 



Equivalently, one can write Eq. OH as: 
(b„ ^ + C„) ij + Ci2;r = /I [Ell (f/ - *) + Ei2;e] 

{B2lf^+C2l)U+ (B22^+C22);r = 
A [ (D21 1 + E21) {U-^')+ (D22^ + E22)^] 
C3lf/ + C324'+ (A33^+B33^+C33)* = 



(36) 



Each sub-matrices can be expressed in terms of the two follow- 
ing functionals: 



C(G) 



= 27r r 

Jo 



G{0 0)Y^'(0)Y^;(0) sin ed6 



J";jG) = 27T 



dV" 
G{(,e)Y"/i9)-^{9)sm0de 



(37) 



(38) 



where Y'"(0) - Y'"(0, 0)e"""* is a normalized Legendre polyno- 
mial. 

In Appendix B, all the coefficient of the sub-matrices are 
made explicit in terms of the function r(^, 0) and its first and 
second order derivatives as well as in terms of the enthalpy of 
the equilibrium model, its first and second order derivatives. 

In the following, we consider ffie Cowling approximation 
thus neglecting the gravitational potential perturbation. The 
ODE system J3U is simplified accordingly and in particular 
reduces to the first order 

3. Tests and accuracy 

The formalism and the numerical method presented in the pre- 
vious section have been tested and the accuracy of the fre- 
quency determinations has been estimated. 

3.1. Tests 

A first test of the method has been performed in the case of ax- 
isymmetric ellipsoids of uniform density. We choose this con- 
figuration because the eigenvalue problem is fully separable us- 
ing the oblate spheroidal coordinates (^, rj, (p) defined as (x - 
a cosh ^ sin 77 sin (^, y = a cosh ^ sin 77 cos 0, z - asinh^cos;;). 
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where < ^ < +00, 0<77<7retO<(^< 2n. The eigen- 
frequencies obtained with this method were compared with the 
eigenfrequencies computed with our general method, S (0) de- 
scribing an ellipse. We found the same frequencies with a high 
degree of accuracy for arbitrary values of the ellipsoid flatness 
between and 0.5. Moreover, as the flatness goes to zero, the 
frequencies were found to converge towards the values given 
by a first order perturbative analysi s in terms of fla t ness. M ore 
details about this test are given in iLi gnieres et al.l (120011) and 
iLi gnieres & Rieutord.(.2004i) . 

The frequencies of axisymmetric p-modes in a self- 
gravitating uniformly rotating N = 3 polytrope that wifl be 
presented in the following sections have been also tested. As 
shown in the previous section, the method involves lengthy 
analytical calculations of the coefficients of the ODE system 
<31> . Terms involved in the non-rotating case have been tested 
by comparing our result with the p-modes frequencies in a 
non-rotating self-gravitating N = 3 po lytrope published in 
dChristensen-Dalsgaard & Mullanlll994l) . The relative error is 
smaller than 10"^ for the ^ = Oto3,«=ltolO modes. In the 
rotating case, we compared our results with the ones obtained 
by solving the same problem but using a different form of the 
starting equations. This alternative system of equations aims at 
including the Coriolis force; thus the variables and the result- 
ing ODE systems to be solved are different. We verified that 
when the terms involving the Coriolis force are omitted from 
the equations, eigenfrequencies presente d in the next sectio n 
are recovered with a very high precision (Rees e et al.L l2006aS. 
For instance, the maximum relative error on the eigenfrequency 
for all the frequencies computed at Q/Q^ = 0.59 has been 
found of the order of 10"^, for a given set of numerical param- 
eters. This last test gives us strong confidence in the method 
and its implementation. 




Fig. 1. Evolution of the frequency relative error, E(L) - \(d(L) - 
(y(L™'"')|/w(L™'"'), as the spatial resolution in latitude is in- 
creased. Two modes labeled ({ - I, n = 1) and {( - 6, n = 8) 
are considered at the rotation rate Q.ID.K - 0.46 with L""*" = 80 
and A^f = 61. 



3.2. Accuracy 

As pointed out bv IClementI ( ll98lUl998l) . accurate numerical 
solutions of the 2D eigenvalue problem are not easy to obtain. 
It is therefore important to estimate the accuracy of our method. 
In the following we first investigate the influence of the spa- 
tial resolution on the eigenfrequencies and then consider other 
sources of errors. 




Fig. 2. Evolution of the frequency relative error, E(Nr) — 
\co(Nr) - w(A^™'')|/w(A^™'') as the resolution in the radial co- 
ordinate is increased. Two modes labeled (f = 1, n = 1) and 
(t = 6, n = 8) are considered at a rotation rate Q./Q.K = 0.46 
with A?™" = 61 and L = 62. 



A relative spectral error E is defined as the absolute value 
of the relative difference between the frequency computed at 
a given resolution and the frequency obtained at the maxi- 
mum resolution considered. Let us first consider the effects 
of the angular resolution. Fig. [Q displays E(L) = \cl){L) - 
w(L™^'')|/aj(L™^'') as a function of L, the truncation order of 
the spherical harmonic expansion, for two axisymmetric modes 
labeled (/' = 1, « = 1) and {{ - 6, n = 8) whose spatial struc- 
tures are dominated by large and small length scales, respec- 
tively (the labeling of the mode wifl be described in the next 
section). The maximum angular resolution is L™^" = 80, the 
resolution in the pseudo-radial coordinate is fixed to Ar =61 
and the rotation rate is Q.IO.K - 0.46. In the same way, Fig.|2 
illustrates the effects of the pseudo-radial resolution by show- 
ing E(Nr) = \oj{Nr) - w(A,"^'')|/w(A™'') as a function of Nr, 
the truncation order of the Chebyshev polynomial expansion, 
for the same modes and rotation rate. The maximum radial res- 
olution is A™"" =61 and the latitudinal resolution is fixed to 
L - 62. In both figures, we observe that the error first decreases 
and then reaches a plateau which means that a better approxi- 
mation of the eigenfrequency cannot be obtained by increasing 
the spatial resolution. The plateau are significantly higher for 
the {( = 6, n - 8) mode than for the (^ = 1, n = 1) mode. We 
verified that this difference is due to the presence of smaller ra- 
dial length scales (rather than to smaller angular length scales). 

Even when the spatial resolution is sufficient, two other 
sources of numerical errors can indeed limit the accuracy of 
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eigenfrequency determination. First, the component of the ma- 
trix L and M being computed numerically they are determined 
with a certain error. Second, even when this error is reduced 
to round-off errors, the accuracy of the algebraic eigenvalue 
solver, the Arnoldi-Chebyshev algorithm, remains limited. 
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Fig. 3. Histogram of 100 frequencies obtained for 100 different 
values of the initial guess of the Arnoldi-Chebyshev algorithm 
randomly chosen in a small interval around w,„ - 12.0547. The 
standard deviation of this frequency distribution cr = 5.6 x 10"^ 
measures the algorithm relative error on the frequency for this 
particular mode labeled ({ - b, n - 8). The spatial resolution 
is A^r = 61 and L- 62 and the rotation rate is Q./Q.K = 0.46. 



Errors on the matrix component that arise from quadratures 
(see Eqs. (I37> and (I38» can approach round-off errors using 
weighted Gauss-Lobatto quadratures. The other source of er- 
ror in the matrix components comes from the computation of 
equilibrium quantities. Indeed, the accuracy of the enthalpy, its 
first and second derivatives and the surface shape, is at best 
limited by the effect of round-off errors on the convergence of 
the algorithm used to compute the polytropic stellar models. 
The effect of these errors on the eigenfrequencies have been 
investigated and appears to be smaller than the effect of the 
Arnoldi-Chebyshev algorithm itself which is now described. 

As any solver in linear algebra, the Arnoldi-Chebyshev 
algorithm amplifies the round-off error that affect the matrix 
components. Thus, the error on the eigenvalue and the as- 
sociated eigenvector is usually much larger than the round- 
off error of double precision arithmetic. The accuracy of the 
Arnoldi-Chebyshev algorithm has been studied in details by 
IValdettaro et al. (2000) in the context of inertial modes in a 
spherical shell where the matrix component are known ana- 
lytically. Theoretically, it can be estimated by computing the 
spectral portrait of the eigenvalue problem [M]X - A[Q]X, 
which shows how small variations of [M] and [Q] affects the 
determination of each eigenfrequencies. In fact, as the itera- 
tive Arnoldi-Chebyshev algorithm requires an initial guess of 
the eigenfrequency, a practical alternative to measure the ac- 
curacy of a frequency determination is to compute frequencies 
for slightly different values of the initial guess. This has been 



done for a large number (100) of initial guess values randomly 
distributed around the eigenfrequency of the {£ - I, n - I) 
and {1-6, n = 8) modes. The histogram in Fig. |3l shows 
the resulting frequencies distribution around a most probable 
mean eigenfrequency. The width of the histogram determined 
by the standard deviation of the distribution provides a measure 
of the algorithm accuracy. The standard deviation cr is equal to 
5.6 X 10"'' for the (^ = 6, n = 8) mode and to 6.2 x 10""' 
for the ({ = \, n = \) mode. The error thus grows with the 
radial order of the mode, this trend being general in our results 
(as in Valdettaro et al., 2000). Moreover, the width of the his- 
togram does not depend on the amplitude of the initial guess 
perturbation provided it is sufficiently small. 

We also observe that, except for the dependence of the {{ - 
1, n = 1) frequency on the angular resolution, the levels of the 
plateau shown in Figs.^and Hare of the same order as the error 
of the algorithm. It means that, in these cases, the changes in the 
matrix component and size associated with the modification of 
the resolution have a similar effect on the frequency as varying 
the initial guess of the algorithm. However, the convergence of 
the (^ = 1, n = 1) frequency at a 10"''* level, much lower than 
the 6.2 X 10"'° accuracy of the algorithm, shows that it is not 
always true and that the spectral error can underestimate the 
true error. 

Although it is too demanding to compute a global accu- 
racy by repeating the statistical study on the initial guess for all 
eigenfrequencies, the relative accuracy on all tested frequencies 
is always better than 2x 10~^ using double precision arithmetic. 

Note that another potential source of error will be discussed 
below when describing avoided crossings between modes. 



4. Results 

The parameter range of the calculations is first presented. Then, 
we describe the method used to label the eigenmodes, the struc- 
ture of the frequency spectrum, some properties of the eigen- 
functions and the differences with perturbative methods. 



4.1. Parameter range 

Self-gravitating uniformly rotating polytropes of index N - 3 
and specific heat ratio Fi o - 5/3 have been computed for ro- 
tation rates varying from Q = up to O/Qj^ = 0.59. In this 
range, the flatness of the star's surface e - I - Rp/Re increases 
from to 0.15. 

Low frequency axisymmetric p-modes have been computed 
for each polytropic model. We started with the non-rotating 
model and computed the { - (0, ...,7), n = (1, ...,«,„„) ax- 
isymmetric p-modes, the largest radial order depending on the 
degree {: n„,a, = 10 for { = (0, 1), «,„«,- = 9 for ^ = (2,3,4) 
and n,„„v = 8 for { - (5,6,7). All these 71 modes were then 
calculated at higher rotation rates by progressively increasing 
the rotation of the polytropic model. In the next section, we ex- 
plain how we could track and label them from zero rotation to 
Q./Q.K = 0.59. 
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Fig. 4. Evolution of all the computed p-modes frequencies from Q = to Q/Q^ - 0.59. The frequencies have been adimen- 
sionalized by {GMIRiY^^ because we expect that the polar radius Rp does not change much as the rotation of the star increases. 
Non-rotating { - (0, ..., 7), n - {\, ..., n,„a,i) p-modes have been followed by progressively increasing the rotation. This mode 
tracking requires special care when an avoided crossing occurs between two modes of the same equatorial parity. The figure on 
the left shows an overview of the frequency evolution while the two right figures display zooms to illustrate avoided crossings 
between the E - \, n - b and ^ = 5, n - 5 modes and the ( = 0, n = 4 and ^ = 4, n - 3 modes, respectively. Although the two 
"interacting" modes have a mixed character near the closest frequency approach, their original properties are recovered after the 
crossing which enables to unambiguously follow and label the modes. This is illustrated in Fig.|5lby considering the spectra of 
Legendre expansion components of the i = 0, n = 4 and t - 4, n = 3 modes at the rotation rates marked by an arrow. Note that 
in the above figures crossings do occur between equatorially symmetric and anti-symmetric modes. In the global view, there are 
two examples of discontinuous frequency changes due to avoided crossing with modes which frequency is not represented on the 
figure. Actually, the /" = 8, n = (1, 2, 3) modes have been displayed in this view to avoid more discontinuous changes. 



4.2. Mode labeling 

In the absence of rotation, modes are identified and classified 
by the three "quantum" numbers n, {, m characterizing their ra- 
dial, latitudinal and azimuthal structure respectively. Because 
of separability, independent ID eigenvalue problems are solved 
for each couple {{, m) and it is then an easy task to order 
the computed frequencies, the order n additionally indicating 
the number of radial nodes of the mode. By contrast, in the 
presence of rotation, independent 2D eigenvalue problems are 
solved for a given | m \ and a given equatorial parity. The com- 
puted modes are then obtained without a priori information 
about their latitudinal and radial structures. Therefore, an im- 
portant issue is whether it is possible to define a meaningful 



classification of these modes. In this present work, we investi- 
gate the possibility to associate unambiguously each mode with 
a non-rotating mode thus identifying it by the three quantum 
numbers n, £, m of the non-rotating mode. Similarly, ClementI 
(■1986 ) followed some equatorially symmetric acoustic modes 
to high rotation rates but in a limited frequency range and using 
low spatial resolution calculations. 

In practice, instead of backtracking modes towards zero ro- 
tation, we started at zero rotation with a mode we are interested 
in and tried to follow it by progressively increasing the rotation. 
We did manage to track all the { = (0, ...,7), n = (1, ...,nmax) 
axisymmetric p-modes from Q = to Q/Q;f = 0.59, a 
global view of the eigenfrequencies evolution being displayed 
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Fig. 5. Evolution of the Legendre components C{{) = max„JU{£,nr)\ of the ^ = 0, n - A and ( - A-, n = 3 modes during 
the avoided crossing shown in Fig.|3 Prior to (Q/Q/f = 0.38) and after (Q/Q/f - 0.55) the avoided crossing, the spectrum of 
Legendre components peaks at a given degree while, near the closest approach (Q/Qa: = 0.46), the double peaks of the spectra 
show the mixed character of the eigenmodes. 



in Fig. 13 (left panel). As explained below, the main difficulty 
comes from avoided crossings between modes of the same 
equatorial parity. 

Zooms in the w - Q plane displayed in Fig. 0] (right pan- 
els) provides two examples of avoided crossings respectively 
between odd {( - I, n - 6 and i - 5, n - 5) and 
even {{ — 0, n - A and i - A, n - 3) modes. Modes 
tends to cross because their frequencies are not affected in 
the same way by the centrifugal force but, as two eigenstates 
of the same parity cannot be degenerate, an avoided cross- 
ing takes place during which the two eigenfunctions exchange 
their characteristics. This exchange of property is illustrated in 
Fig. 13 in the case of the (^ = 0, n = 4), (t ^ 4, n = 3) 
crossing. A mean Legendre spectrum is displayed before, near 
the closest frequency separation and after the avoided cross- 
ing. The mean Legendre spectrum of a field U is defined as 
C(() - maxn^\U{l,nr)\lmax\U(i,nr)\, where U{l,nr) are the 
components of the Legendre/Chebyshev expansion, n^ being 
the degree of the Chebyshev polynomial. The quantity Citf 
thus represents the largest Chebyshev component for a given 
value of I normalized by the maximum over all spectral com- 
ponents. The mean Legendre spectra peak at one characteris- 
tic degree before and after the avoided crossing, thus showing 
that the modes recover their original properties after the cross- 
ing and therefore can be unambiguously recognized. Up to the 
fastest rotation considered, the ( = 0-1, « = 1 - 10, m = 0, p- 
modes undergo a limited number of avoided crossing and could 
be followed unambiguously. 

It remains that near the crossing the labeling is somewhat 
ambiguous. First, it is difficult to define a criterion to assign 
a label. Here, we mostly use the degree at which the mean 
Legendre spectrum reaches a maximum. But it occurred that 
the two interacting modes peaks at the same degree in which 
case we determined the location of the smallest frequency sep- 
aration. Second, as shown by Fig.|5] a more fundamental prob- 



lem is that a single label can not reflect the mixed nature of the 
eigenfunction. 

Another issue related to avoided crossings concerns their 
influence on the accuracy of the eigenfunction computation. 
Indeed, if a large-scale well resolved eigenfunction undergoes 
an avoided crossing with a small scale unresolved mode, the 
accuracy of the eigenfunction determination will be affected. 
The effect on the frequency accuracy should be small as the 
frequency gap induced by the avoided crossing of two modes 
of well separated length scales is small. But, right at the closest 
approach, the eigenfunctions will be much affected. At zero ro- 
tation, the highest degree mode present in our frequency range 
is ( = 5l,n = I. Thus, if one of the low degree mode that we 
computed undergoes an avoided crossing with a mode of such 
a high degree, the high degree mode should be resolved to en- 
sure an accurate determination of the eigenfunction of the low 
degree mode. 

4.3. The Structure of the frequency spectrum 

The effect of the centrifugal force on the acoustic frequency 
spectrum of axisymmetric modes is investigated. The mean 
modifications of the spectrum are first commented. Then, we 
investigate how regularities in the frequency spacings evolve 
with rotation. Finally, differences between equatorially sym- 
metric and anti-symmetric modes are outlined. 

4.3.1. Global spectrum evolution 

Figure|6lcompares the frequency spectrum of the ^ = - 7, n = 
1 - n„ax, in = modes at Q = (upper panel) and at Q./Q.fc = 
0.59 (lower panel), the height of the vertical bars corresponding 
to the degree £ of the mode. It appears that the centrifugal force 
induces a mean contraction of the frequency spectrum. This is 
expected as the decrease of the sound speed and the increase 
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Fig. 6. Frequency spectrum of ^ = 0-7, n= 1- nmax, m = modes at Q = (top panel) and Q./Q.K = 0.59 (bottom panel). The 
degree number {! associated with the frequency is shown by the height of the vertical bar. 



of the star volume induced by the centrifugal force both tend to 
lessen the frequency of acoustic modes. 

To illustrate this effect, one could think of a spherically 
symmetric decrease of the effective gravity. If one considers a 
homologous series of spherical models of increased volume V, 
the decreasing rate of the frequencies Alj/o) is -{l/2)(AV/V), 
as the normalized frequencies (L)/{GM/R^y^^ remain constant. 
For non-homologous spherically symmetric changes, Acj/o) is 
asymptotically equal to -A(ln J '' dr/c) for high order modes 
verifying the following asym ptotic formula valid for low de- 
gree and high order p-modes (lTassoulLll980l) : 



CO = 



Jo 



-(n + {e+l/2)/2 + a) 



(39) 



where 1 / J ^ is the sound travel time along a stellar radius 
and a is a constant. When, as in these two previous cases, 
Aw/tLt does not depend on the frequency, the concentration of 
the frequency spectrum is homothetic. 

This is clearly not the case here since the frequencies cross 
each other (see Fig.0}. But there is still an average contraction 
rate which is of the order of -(1 /2)(Ay/y), where now V is the 
volume of the centrifugally distorted star. In addition, the con- 
traction rates of individual frequencies appears to be comprised 



between the logarithmic derivative of the sound travel times 
computed respectively along the polar and equatorial radii: 



5n|ln r"^]<-5n(lno.)<5n|ln 



Jo c 



(40) 



Another interesting property is that, at small rotation rates say 
Q./D.K < 0.05, the contraction rate 5n(lnw) tends to be inde- 
pendent of { and n for the large degree modes H > 3. This sug- 
gests that an asymptotic regime exists for modes with horizon- 
tal wavelengths smaller than the dominant length scales of the 
centrifugal distortion. In this regime, the contraction rate has 
a constant value which is not equal to -(l/2)(Ay/y) and that 
would be interesting to determine. We alread y found such be- 
havio ur in the case of homogeneous ellipsoids ( iLignieres et all 
I2OOII) where a perturbative analysis shows that the contraction 
rate of axisymmetric modes is constant for high { and n and 
that it can be related to the increase of the ellipse perimeter. 

Nevertheless, for the low degree modes £ < 2 below 
Q.IQ.K ~ 0.05 and for all modes at higher rotation rates, 
5n(lnw) depends on { and n. This differential effect modi- 
fies the structure of the frequency spectrum as the rotation in- 
creases. 
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Fig. 7. Regularities in the frequency spacings of axisymmetric (m - 0) modes. The large frequency separation between modes 
of consecutive order A„ = a>„( - (jl>„-ic, the frequency separation between t + 2 and I modes, A2/ = w„/+2 - w„/, and the small 
frequency separation 6 - IS.,, - ^2/ are displayed as a function of the radial order n for four different rotation rates, (a) = 0, 
(b) Q.IQ.K - 0.32, (c) Q./D.K - 0.46 and (d) Q./Q,k = 0.59. We plotted the opposite of the small frequency separation -6 for the 
clarity of the figure. Continuous lines have been drawn between frequencies of the same degree {. 



4.3.2. Regular frequency spacings 

In a non-rotating star, the frequency spectrum presents some 
regular frequency spacings which can be accounted for by an 
asymptotic theory in the high frequency limit w — > 00. The 
asymptotic formula ( I39> . valid for low degree and high or- 
der modes, shows that the large frequency separation between 
modes of consecutive order n , A„ = (^n+i,e - <^n,t, does not 
depend on ( and n and is equal to ttj J ^. A more detailed 
asymptotic analysis also shows how the so-called small fre- 
quency separation 6 - con+\,t - w„/+2 vanishes as a function 
of the frequency. Although our calculations are restricted to the 
low frequency part of the acoustic spectrum, we do observe a 
clear tendency towards these asymptotic behaviors in the non- 
rotating case. We can therefore investigate whether these prop- 
erties are modified by rotation. 

Figure0presents the large frequency separation A„ and the 
frequency separation between consecutive modes of the same 
order and parity: 



A2/ 



(^n,(+2 - W„/, 



(41) 



as a function of the radial order n for four different rotation 
rates, (a) Q = 0, (b) Q./D.K = 0.32, (c) Q/Q^^ = 0.46 and 
(d) Q./Q.fc = 0.59. As in the previous figures, the frequencies 
are adimensionalized by {GM/R^^y^^. Continuous lines have 
been drawn between frequencies of the same degree L We first 
observe that the large frequency separation tends to be inde- 
pendent of n and { at all rotation rates. In accordance with the 
mean contraction of the frequency spectrum mentioned above, 
the large frequency separation decreases with rotation. It is al- 

pR rR 

ways comprised between jil L '' dr/c and n/ L ' dr/c. 

The dispersion of the large frequency separations around 
their mean value also has an interesting evolution with rotation. 
In the non-rotating case, the dispersion reflects a regular depar- 
ture from the asymptotic limit. It is larger for high degrees and 
monotonically decreases with frequency (see Fig.^). In the ro- 
tating cases, the dispersion is not as regular The largest depar- 
tures, some of which are most clearly visible in Fig.Et, can be 
attributed to an ongoing avoided crossing. The residual disper- 
sion is iiTegular and decreases with rotation. At Q/Qa: = 0.59, 
if we exclude all « < 4 values from our sample, the mean large 
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frequency separation < A„ > is equal to 1.095{GM/R^y^^ and 
its standard deviation is 0.017 < A„ >. 

We now consider the small frequency separation 6 = 
A„ - Aa.f. As expected, in the absence of rotation the small 
frequency separation tends to vanish as n increases. But, at 
n/Q.K = 0.32, the small frequency separation no longer de- 
creases with n for some values of { and for the higher rotation 
rates it becomes nearly constant. At the same time, the A2/ sep- 
aration becomes more and more uniform as rotation increases. 
As shown in Fig. Eb, A2/ becomes approximatively constant 
with n first for low degree modes while it still increases with 
frequency for high degree modes. In addition, equatorially an- 
tisymmetric modes reach this new regime at a lower rotation 
rate than the equatorially symmetric modes of similar degree. 
This is illustrated in Fig. 01 by the A2/=4 curve which still re- 
mains above the mean A2,f value while the A2.f=5 separation 
already collapses with the other curves. At Q./D.K - 0.59, if we 
exclude all n < 4 values from our sample, the mean frequency 
separation A2,f is equal to 0.381 (GM/R^y^- and its standard 
deviation is 0.033 < A„ >. 

As a consequence of the near uniformity of A„ and A2,f, 
the frequencies of low degree and high order modes can be 
approximated by the following expressions: 



C^nC 



n6n + p5( + a^ if { - 2p 
n6n + pSf + a^ if ( - 2p + \ 



(42) 



where 6„ -< A„ >, 6c -< A2/ >, c^ and a" only depend on the 
equilibrium model. Using a reference frequency to determine 
the a constants (the i - Q, n - E frequency for a^ and the 
£ - \, n = 8 frequency for a ), we computed the root mean 
square error -\J1/NYj(<^ - w)^ and the maximum error made in 
using the approximate expressions (I42t . For a frequency subset 
containing the n > 4 and { < 5 modes, the rms error is 0.0116,, 
while the maximum error amounts to 0.055,,. Both errors are a 
very small fraction of the large separation which shows that Eq. 
i42> yields useful approximations of the frequency spectrum. 

4.3.3. Equatorially symmetric versus anti-symmetric 
frequency spectra 

We have seen that the regular frequency spacings A„ and A2/ 
have similar values for symmetric and anti-symmetric modes 
with respect to the equator. The evolution of the equatorially 
symmetric and anti-symmetric frequency spectra are neverthe- 
less quite different. Indeed, considering two modes of similar 
frequency but of opposite equatorial parity, the frequency of the 
symmetric mode generally decreases faster with rotation than 
the frequency of the antisymmetric modes. The consequence 
is that the frequency separation between modes of consecutive 
degree (and thus of opposite parity) A; - (d„(+] - (l)„[ tends to 
increase when £ is even and to decrease when £ is odd. The fre- 
quency separation A; can even become negative which implies 
that, contrary to the non-rotating case, frequencies of a given 
order n do not increase monotonically with the degree £. This 
striking modification of the usual frequency ordering is appar- 
ent on Fig. |S] where the (£ = 2, n) frequencies are smaller than 
the {£ - 1 , n) frequencies for all the order n that we calculated. 



that is n = (1, ..., 10). In the same way, the (£ - 4, n) frequen- 
cies are smaller than the {£ = 3,n) frequencies if n > 3, and 
again the {£ = 6, n) frequencies are smaller than the {£ = 5, n) 
frequencies if « > 5. 

4.4. Equatorial concentration 

In this section, we shall focus on the most notable effect of 
the centrifugal force on the eigenmodes, namely the equato- 
rial concentration and consider its consequences on the mode 
visibility. Note that Clement (1981) also reported an equato- 
rial concentration of the equatorially symmetric modes that he 
calculated. 

Figure [S] shows this effect on the {£ - A,n - A) mode. 
Contours of the amplitude of the Lagrangian pressure pertur- 
bation are plotted in a meridional plane for increased rotation 
rates, (a) Q = 0, (b) D.IQ.K = 0.32, (c) D.ID.K = 0.46 and (d) 
Q./D.K = 0.59. We observe that the number of nodes increases 
along the equatorial radius and decreases along the polar one. 
Along the surface, the number of nodes remains equal to £ be- 
fore Q/Q;f = 0.59 where additional nodes appear The equato- 
rial concentration is clearly seen in the outermost layers. 

In Fig. |5] the equatorial concentration is shown for other 
modes including the lowest and highest degree modes of our 
sample as well as symmetric and anti-symmetric modes. The 
latitudinal variation of the mode amplitude is displayed at 
the surface for the following modes, (a.) £ - 0, n - 1, (b) 
/>= 1,„= l,(c)/' = 6,n = 8, (d)^ = 7,« = 8. In each case, the 
equatorial concentration grows with rotation. At the largest ro- 
tation rate, symmetric modes are maximal at the equator while 
anti-symmetric modes peak at small latitudes since they must 
vanish at the equator The contrast between these maxima and 
the polar amplitude is strong. 

The equatorial concentration reveals a modification of the 
resonant cavity of the acoustic waves. In particular, the reduc- 
tion of the volume of the resonant cavity should tend to increase 
the frequency. The equatorial concentration seems also to be 
associated with the near-uniformity of the frequency separation 
A2,f . At small rotation rates, the concentration is not completed 
and A2/ is clearly not constant. Whereas, at the largest rota- 
tion rate, all modes are concentrated near the equator and A2/ 
is nearly uniform. 

Besides its effect on the frequency spectrum, the equatorial 
concentration of eigenmodes should induce a profound modi- 
fication of the mode visibility as compared to the non-rotating 
case. The photometric mode visibility is determined by the in- 
tegration over the visible part of the star's perturbed surface of 
the radiation intensity perturbations associated with a particular 
pulsation mode. Rigorous calculations of photometric visibili- 
ties are beyond the scope of the present paper as they require 
non-adiabatic calculations of the oscillation modes and stellar 
atmosphere models (e.g. Daszyiiska-Daszkiewicz et al. 2002 ). 
But we can still determine the effects of averaging the pertur- 
bations over the visible surface which have a direct impact on 
the visibility. The disk-averaging factor is defined as: 



D(i) 



nRl6To iis,^ 



6T{6, 0)dS ■ Ci 



(43) 
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Fig. 8. Isocontours of the £ = 4,n - 4 mode amplitude in a meridional plane as a function of the rotation (a) Q = 0, (b) 
Q./Q.K = 0.32, (c) Q./Q.K - 0.46 and (d) Q.IQ.K - 0.59. The amplitude is normalized with the maximum of its absolute value. 
Continuous lines correspond to positive amplitudes, dashed lines to the zero amplitude and dotted lines to negative amplitudes. 
At zero rotation, the angular distribution is given by the Y-liff) Legendre polynomial while the radial distribution is characterized 
by the surface concentration of the amplitude and the presence of n = 4 nodes in the inner part. For larger rotation rates, the 
largest amplitudes concentrates toward the equator. We also note that the number of radial nodes decreases along the polar radius 
while it increases along the equatorial radius. 



where i is the inclination angle between the line-of-sight and 
the rotation axis, Cj is a unit vector in the observer's direction 
and 6T is the spatial part of the Lagrangian temperature per- 
turbation at the stellar surface, 6T being proportional to the 
velocity divergence x in the approximation of adiabatic per- 



turbations. The mode amplitude is normalized by 6T() the root 
mean square of the perturbation over the whole stellar surface 



^ro 



-{K 



ST^iO, (l>)dS 
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(44) 
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Fig. 9. The mode amplitude at the surface of the polytropic model as a function of the rotation rate (&) £ - 0,n - 1, (b) 
{ = \,n - 1, (c) ^ = 6, n = 8, (d) ^ = 7, n = 8. The amplitude is normalized with the maximum of its absolute value. While the 
angular distribution is given by the corresponding Legendre polynomial T^{9) in the absence of rotation, the oscillation amplitude 
progressively concentrates towards the equator (9 - n/T) as rotation increases. 



and the visible surface S ,, has been normalized by nR^, the vis- 
ible surface of a star seen pole-on. With these normalizations 
the disk-averaging factor of a radial mode seen pole-on is unity. 

In the absence of rotation, the surface distribution of modes 
is determined by a unique spherical harmonic and the disk- 
averaging factor takes a simple analytical form (Dziembowski, 
119771) . For even degree and for £ - 1, the disk-averaging fac- 
tor varies with the inclination angle as the Legendre polyno- 
mial Y"\i) while it vanishes altogether for odd degree £ > 3. 
For rotating stars, the method of the calculation is detailed in 
Appendix C. Note that for modes which are equatorially anti- 
symmetric and axisymmetric, the disk-averaging factor has 
also a simple dependency on the inclination angle as it is pro- 
portional to cos(/). 

Figure ^| shows the disk- averaging factor of various ax- 
isymmetric modes as a function of the inclination angle. The 
non-rotating case is displayed in Fig.llOh where /" = to f = 7 
modes are considered. We recall that, at Q = 0, modes of 
different radial orders but same £ and m have the same sur- 
face distribution. This is not true in the rotating case and, at 
Q./Q.K - 0.59, Fig.llOb andllOt present the disk-averaging fac- 



tor for modes of the same degree numbers but for two different 
radial orders n - I and n = 8, respectively. Note also that the 
disk-averaging factor was allowed to take negative value for 
the clarity of the figure although it is its absolute value which 
is relevant for the mode's visibility. Figure [TO] shows that ro- 
tation strongly modifies the dependency of the disk- averaging 
factor on the inclination angle as well as on the degree number. 

Fig. not shows that for all n = 8 equatorially symmetric 
modes the absolute value of disk-averaging factor tends to in- 
crease with the inclination angle. This is due to the equatorial 
concentration of these modes (see for example the surface dis- 
tribution of the (^ = 6, « = 8) mode at Q/Q.^ - 0.59 shown in 
Fig-Et)- This tendency is less pronounced for the n = 1 sym- 
metric modes shown in Fig. llOb (except for the { - 2 mode) 
although these modes are also equatorially concentrated. This 
is due to a cancellation effect between positive and negative 
perturbations concentrated near the equator as illustrated by the 
surface distribution of the {£ - Q,n - 1) mode in Fig.|9^. The 
non-rotating case strongly differs since the absolute value of 
the disk-averaging factor for even degree £ > Q modes does not 
vary monotonically with the inclination. Indeed, they have £12 
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Fig. 10. The disk-averaging factor D(i) is shown as a function of the inclination angle / for various axisymmetric modes at two 
different rotation rates O = (a) and Q./Q.K = 0.59. The degree of the modes varies from ^ = to ^ = 7. In the rotating case, 
the surface distribution also depends on the order of the mode. Two values n = 1 (b) and n = 8 (c) have been considered at 
Q/Q^ = 0.59. 



(a) .= 





Fig. 11. The variation of the disk-averaging factor as a function of the degree { is shown at two fixed values of the inclination 
angle (/ = and / = 7r/2). The three curves on each figure correspond to the non-rotating case and, at Q./Q.K = 0.59, to two 
different orders, « = 1 and « = 8. The absolute value of disk-averaging factor has been rescaled by its maximum value among 
the different degree considered to outline its dependency on the degree number. In sharp contrast with the non-rotating case, the 
disk-averaging factor at Q/Q/f = 0.59 does not show a strong decrease with {. 



nodes between and n/2. For odd { modes, the disk-averaging surface distribution onto the Legendre polynomial FJ* is not 

factor is also modified by rotation since it no longer vanishes zero for £ > 3 modes. 

for e > 3. This occurs because the projected elementary sur- :„ non-rotating stars, the cancellation effect between pos- 

faces dS ■ ei are no longer symmetric with respect to the ob- itive and negative perturbations results in a rapid decrease of 

server's direction and because the projection of the eigenmode the disk-averaging factor as the degree { of the mode increases. 
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Consequently, modes above a certain degree E > 3 - 4 are 
not expected to be detectable with photometry and are there- 
fore not included when trying to identify the observed frequen- 
cies. As shown in Fig. ^2 this property must be reconsidered 
for rapidly rotating stars. The absolute disk-averaging factor 
normalized by its maximum value over the degree considered 
Q < { < 1, \D(i)\/max(\D{i)\, is plotted as a function of d for 
two fixed values of the inclination angle, / = on Fig.lllh and 
/ - jt/2 on Fig.lllb. The three curves correspond to Q = and 
to Q,/Q.K - 0.59 for the n = 1 and n = 8 modes, respectively. 
In sharp contrast with the non-rotating case, the disk-averaging 
factor has no tendency to decrease above ( - 2. Again, this 
can be explained by the equatorial concentration as modes of 
different degree have a similar surface distribution. 



4.5. Comparison with perturbative metliods 

According to the perturbative anal y sis, ce ntrifugal eff'ects ap- 
pear at second order in D. (Saio, 1981). To determine the 
second-order perturbative coefficient from our complete calcu- 
lations, we performed a serie of calculations for small rotation 
rates (Q = 0, 1.8 x 10"^ 1.8 x 10-^,4.6 x 10^^0.09 x 10"-, ... 
in units of Q.k)- From them, we determined the second-order 
perturbative coefficient, denoted wi, as the limit of the ratio 
(w(Q) - a>())/ii^, where wq denotes a non-rotating eigenfre- 
quency. Thus the approximate frequencies valid up to the sec- 
ond order in Q. read Wpgrt = (^o + cj\^^, where the frequencies 
are in units of {GMIE?y^^ and the rotation is in unit of Q^:. To 
assess the range of validity of the second order perturbative ap- 
proach, we compared these approximate frequencies to the "ex- 
act" frequencies. In Fig.[2J the relative differences between the 
two calculations, (w - copQ^)lco, is plotted as a function of the 
rotation rate for the £ -Q-2, n - 1- 10 modes. The departures 
computed for the other modes, { - (3, ...,7), n - (1, ...,«,„„), 
are smaller than the extremal differences shown in Fig. 1121 and 
are not displayed for the clarity of the figure. The relative differ- 
ences are generally larger for low degree modes and, for small 
rotation rates, are a monotonic function of the radial order n (an 
increasing function for the £ - Q -2 modes shown in Fig. ll2> . 
As mentioned before the low degree modes seem to be sensi- 
tive to the precise form of the distortion which occurs at sim- 
ilar lengthscales. As rotation increases, it appears that higher 
than second order effects are important to describe the effect 
of the centrifugal distortion on these modes. The second order 
approximation is much better for large ( modes which are sen- 
sitive to global distortion properties. 

As compared to the observational uncertainties on the fre- 
quency determinations, the error made in using second order 
perturbative methods becomes rapidly significant as rotation 
increases. For a ratio Q./Q.K - 0.24, corresponding to a typical 
6 Scuti star with an equatorial velocity of 100 kms"', the max- 
imum absolute difference amounts to lljuHz which is much 
larger than typical observational uncertainties. 

The fact that in our frequency sample the absolute differ- 
ence increases with frequency suggests that departure from the 
perturbative approach could be detectable for moderately ro- 
tating stars pulsating in high order modes. In our data limited 
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Fig. 12. The relative difference between exact frequencies and 
their second order perturbative approximation (second order in 
terms of the small parameter fi/Q^:) namely, doj/co where 6a) - 
a> - Wpert is displayed as a function of the rotation rate for the 
^ = - 2, n = 1 - 10, m = modes. 

to n < 10, the largest relative difference is 2.87 x lO""' for a 
Q./Q.K = 0.139 rotation corresponding to a solar-type star with 
an equatorial velocity of 60 kms"' . If we assume that the same 
relative difference holds for high order p-modes in the range of 
2000//Hz generated in convective envelopes of these stars, the 
absolute difference is S.lfjHz for a typical 2000/iHz frequency. 
This would be also easily detect able given the observational 
uncertainties rBouchv et al., 2005). However, a firm conclusion 
should await a direct comparison between the perturbative ap- 
proach and the complete calculation for p-modes generated in 
the convective envelope of rotating solar- type stars. 

5. Discussion and conclusion 

A new non perturbative method to compute accurate oscillation 
modes in rotating stars has been presented. The accuracy of the 
computed frequencies has been obtained by testing the effect 
of the different parameters of the numerical method. Then, the 
effects of the centrifugal force on low frequency axisymmetric 
acoustic oscillation modes have been investigated in uniformly 
rotating polytropic models of stars. Seventy-one low degree 
£ < 7 and low order « < 10 modes have been first determined 
at zero rotation and then tracked at higher rotation rates up to 
Q./£Ik = 0.59. 

In the frequency and rotation ranges considered in this pa- 
per, the zero rotation quantum numbers £ and n have been used 
to label the modes. This labeling turned out to be meaningful 
since we found regular frequency spacings between modes of 
the same degree and consecutive orders, A„, and, within the 
subsets of modes of the same equatorial parity, between modes 
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of the same order and consecutive degree, A2/. We noted how- 
ever that near avoided crossings, when the eigenfunction is a 
mix of the two "interacting" modes, a unique label cannot re- 
flect the actual eigenfunction structure. Although successful in 
the frequency and rotation ranges considered, it remains to be 
proved that this labeling can be performed in practice at higher 
rotation and at higher frequencies. Indeed, the main difficulty 
of the labeling procedure arises from the avoided crossing be- 
tween modes of the same equatorial parity and such cross- 
ings will be more frequent as the eigenfrequency density in- 
creases with the frequency. The coupling between modes is 
also stronger at higher rotation rates. It might then be neces- 
sary to investigate other tools than the mean Legendre power 
spectrum to characterize the modes. 

The study of the frequency spectrum showed a quite unex- 
pected result, namely that, at the highest rotation rates, a new 
form of organization sets in after the zero rotation asymptotic 
structure of the spectrum has been destroyed. In the absence of 
rotation, the asymptotic theory is directly related to the spher- 
ical symmetry of the stars and ultimately to the integrability 
of underlying ray dynamics. In the presence of rotation, the 
eigenvalue problem is not fully separable and the underlying 
acoustic ray dynamics is most probably not integrable. Then, 
the regular spacings observed at high rotation rates were not 
really expected. They might be the sign of a near-integrable ray 
dynamic rather than a chaotic system. These aspects will be in- 
vestigated through a ray dynamic study of rotating polytropic 
models of stars. 

Most importantly for asteroseismology, the existence of 
regular spacings in the spectrum can potentially provide tools 
for the mode identification in rapidly rotating pulsating stars. A 
complete acoustic frequency spectrum including m + modes 
and the effects of the Coriolis acceleration should however be 
computed and analyzed to assess the practical relevance of 
these regular spacings. 

Apparently, there is a relation between the new spectrum 
structure and the equatorial concentration of the mode ampli- 
tudes. A consequence would be that this spectrum structure 
does not apply to the whole spectrum. Indeed, sufficiently high 
degree modes should still be of the whispering gallery type 
(e.g. Rieutord 2001). Then, being so different from equatori- 
ally concentrated modes, they are not expected to follow the 
same regular spacings. 

Another interesting issue is the difference between the 
modes of different equatorial symmetry. We have seen that 
although the structure of the symmetric and anti-symmetric 
spectra is similar, the frequency spectrum of the symmetric 
modes as a whole seems to evolve independently from the anti- 
symmetric spectrum. The equatorial symmetry also influences 
the "strength" of the avoided crossings measured by the fre- 
quency separation at the closest frequency approach. As illus- 
trated in Fig.|3](right panels), avoided crossings between sym- 
metric modes are always stronger than avoided crossings be- 
tween anti-symmetric modes since they remain further apart. 

Modes undergoing an avoided crossing are particular be- 
cause they have close frequencies and similar eigenfunctions. 
As a consequence, both can be excited to observable levels by 
some excitation mechanism. They are therefore good candi- 



dates to explain the occurrence of close frequencies in observed 
spectrum (Bre2er& Pamvatnvkh, 2006) as well as the associ- 
ated amplitude variations induced by beating between the two 
close frequencies. 

The most striking effect of the centrifugal force on the 
eigenfunction is the equatorial concentration of the mode am- 
plitude. Again, the study of the ray dynamics should help spec- 
ify the conditions in which the sound waves stay focused in the 
equatorial region. As compared to the non-rotating case, the 
equatorial concentration strongly modifies the integrated light 
visibility and in particular its variation with respect to the mode 
degree and the inclination angle. Interestingly enough, our re- 
sults showing a global increase of the disk-integration factor as 
the star is seen equator-on are compatible with observations of 
6 Scuti pulsations which also suggest an inc rease of the pul- 
sation amplitudes with / JSuarez et al.[ 120021) . Another finding 
of practical interest is that, for rapidly rotating stars, the can- 
cellation effect of the disk averaging does no longer sharply 
decrease with the degree of the mode and also varies with the 
order of the mode. Realistic calculations of the mode visibility 
including non-adiabatic calculations of the oscillation modes, 
stellar atmosphere models as well as the gravity and limb dark- 
ening effect will however be needed to draw firm observational 
conclusions. 

Finally it should be reminded that the omission of the 
Coriolis force did not allow a complete treatment of the ro- 
tational effects. However, the effect of the Coriolis force van- 
ishes for sufficiently large frequency (as the time scale of the 
Coriolis acceleration 1 /O becomes much larger than the pulsa- 
tion period) while the modification of the equilibrium model by 
centrifugal force affects all frequencies. Therefore, the results 
presented here should be at least useful for the high frequency 
part of the acoustic spectrum in rotating stars. In a companion 
paper (Reese et al. 2006b), we extend the present results by tak- 
ing into account the Coriolis acceleration which, among other 
things, allows us to specify the domain of validity of perturba- 
tive calculations. 
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ViE' + V2E^ + ViE^ is 
written in the natural basis £, = dOM/dx' or the conjugated 
basis E' verifying Ej ■ E-i - 6ij, g'" are the components of the 
metric tensor and | § | is the absolute value of metric tensor 
determinant. 

From these expressions, we derived the form of the following 
operators: 



go ■ V = (glEi + glE2) ■ di E' = gld^ + gldg 
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r^V ■ ( . A«) = {r'Al,)d^ + (r'A'a)de + r^V ■ Ao[id] (A.6) 



where Ag^ represent the horizontal part of the Laplacian in 
spherical coordinates: 
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We recall that: 

r = (l-eK+A(0 {S(0)-l+e) 



(A.7) 



(A.8) 
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(A. 10) 
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where S(6) describes the stellar surface. 



Appendix A: Linear operators in spheroidal 
coordinates 

Let us now express the linear operators involved in Eqs. J22> . 
( I23> . (I24> . using the spheroidal coordinates given by Eq. ( I25> . 
We need the general expression of the divergence 



Appendix B: Coupling matrix 

The components of the sub-matrices which define the ODE 
system ( I36t are specified below using the functionals I™^, and 
J™, defined in Eqs. OtJ and (|38}: 



A33 
Bii 






(B.l) 
(B.2) 



v-v^^^di {^f[7]v') 



(A.l) B21 l'^f,{h)-i7f, \2h2+h,^^ 

80J 
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Cii 
C12 

C21 



C, \h 



— ; ' ^n 



W «i 



C(/Z4)-C(2/Z2) 
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hid( 



(B.4) 



The components of the ODE, given by Eqs. (IB.U to JB.18> . 
can then be expressed in terms of the enthalpy and its deriva- 
tives, H^, Hg, H^g, H^^Hgg. This has been done in order to min- 
(B.5) imize the numerical error in the calculation of these compo- 
nents. The most useful expressions are: 

(B.6) 
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where { - m + 2k, {' — m + 2k' when applied to the S,], vector 
and ^ = m -H 2/fe -H 1, r = m -H 2A:' -H 1 for H,;,. / ^ 

For polytropic model of index A^, the quantities describing c'f | ~ 
the equilibrium can be expressed in terms of the dimensionless 
enthalpy H as follows: 



^0^ 



r 

T [/!l//ff - I^Hi^g + (d(hi 

-2hi/h3)H^+ [2h2lh-d^h2)He] 
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where A is such that 
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^-^Ihj (H^gH^ - H^^Hg) + (h2d(hi 

-hid(h2) H^ - d^hiH^Hg + d^h2Hl\ (B.25) 
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where he andpc are the dimensional enthalpy and density at the q j^^ - — —h2 

center of the polytropic model. ''i ''i 



(B.32) 
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we deduce that 



Appendix C: Calculation of the disk-integration 
factor 

According to the definition of the disk-integration factor, Eq. 
( I43> . we are led to calculate integrals of the following form: 



jf G(e,4>,i)dfidcf> 



(C.l) 
(C.2) 



G = rA{e, 4>, OW(0)e'""^ 
where 

A(9, 4>, i) = rr^E^ ■ a = Jr^+rje' ■ a 

= r (sin0cos0sin/ + cosflcosi) ^■ 
rg (sin 6 cos / - cos cos sin /) 



(C.l 3) 



(C.14) 



(C.15) 

d d 

cos; — (rsin0) - sin/cos0 — (rcosO) (C.16) 
du do 



where /i = cos(0) and F{0, <p) is the surface distribution of the 

eigenfunction obtained in the coordinate system ( I25> in which , s j » »i, -^ . j- i ,. tu c 

° -' ' — ' where e denotes the unit vector perpendicular to the surface, r 

the polar axis is the rotation axis. The integral is most simply . i i . j . .u . c >- i t-. ». j 

^ , . . .° , ^ Z andrg are calculated at the star surface (, - I. Thus the depen- 



dency of G on /, (p and 6 can be specified as follows: 



calculated in the coordinate system in which the polar axis is 
aligned with the direction of the observer. This coordinate sys- 
tem results from a rotation of angle / around the y axis of the G - A{6) cos /e'""^ - B{6) sin / cos (pe 
original coordinate system, the new angular variables being de- 
noted 6' and 0'. To express G in these coordinates, we use the where 

formula relating the spherical harmonics in both systems: d 

A ^ r— (r sin 0)^^(0) 
+f do 

Yf(9,4>)^ Z dUi)Yr(0\^') (C.3) 5 ^ ^d_ (^^^^,^^(,) 

m'=-t do 

where ^^„„,(0 do not generally have a simple form jEdmondsl It follows that 

Il960l) . Then, using the spherical harmonic expansion of G, we 

obtain: 
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(={) m=-tm'=-t 



Y,Gl(i)Yl'(e,cp) 
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(C.4) 
(C.5) 



so that the integral now reads: 

I/2n = /,„_! + I,„ + I,„+i where 



'm— 



Then, integrating over the longitude 0', from to 2n, the terms 
involving Y'"' {9' , (p') vanish if m' ^ 0. It follows that 



/,„ = cos/A„,(/) 
sin / 

sin; 



+0O +e 

f={) m=-t 

where we used the following relations. 
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(C 8) ^^^ "' terms being defined accordingly. 

Note that for modes which are equatorially anti-symmetric 
and axisymmetric (m - 0), Ao(0 = JoA^Y^(i) and Ai{i) = 
A_i(0 = 0, thus the integral / reduces to: 



(0 
1 



if i is even and £ ^ 0, 
if^ = 0. 



(-1)^TT77ITT if^isoddand^^ 1, 



(C.9) 



(CIO) 



/ = 4;r -0rAQ cos(0- 



(C.27) 



2.4.. .(€+1) 



.2 iff=l. 

Because of axial axial symmetry, the function to integrate reads 
F(O,4>)^W(0)e'""''. (C.ll) 

Then, from the expression of the vector dS at the star's surface: 
dS = dgOMxd^OMdOdip = EixE^dedcp = ylgE^dOdcpiC.U) 



